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Magnetoelectric responses are a fundamental characteristic of materials that break time-reversal 
and inversion symmetries (notably multiferroics) and, remarkably, of "topological insulators" in 
which those symmetries are unbroken. Previous work has shown how to compute spin and lattice 
contributions to the magnetoelectric tensor. Here we solve the problem of orbital contributions 
by computing the frozen-lattice electronic polarization induced by a magnetic field. One part of 
this response (the "Chern-Simons term") can appear even in time-reversal-symmetric materials 
and has been previously shown to be quantized in topological insulators. In general materials 
there are additional orbital contributions to all parts of the magnetoelectric tensor; these vanish in 
topological insulators by symmetry and also vanish in several simplified models without time-reversal 
and inversion whose magnetoelectric couplings were studied before. We give two derivations of the 
response formula, one based on a uniform magnetic field and one based on extrapolation of a long- 
wavelength magnetic field, and discuss some of the consequences of this formula. 

PACS numbers: 73.43.-f, 85.75.-d, 73.20.At, 03.65.Vf, 75.80.+q 



I. INTRODUCTION 

Understanding the response of a solid to applied mag- 
netic or electric fields is of both fundamental and applied 
interest. Two standard examples are that metals can 
be distinguished from insulators by their screening of an 
applied electric field, and superconductors from metals 
by their exclusion of magnetic field (the Meissncr effect). 
Magnetoelectric response in insulators has been studied 
for many years and is currently undergoing a renaissance 
driven by the availability of new materials. The linear 
response of this type is the magnetoelectric polarizabil- 
ity: in "multiferroic" materials that break parity and 
time-reversal symmetries, an applied electric field cre- 
ates a magnetic dipole moment and a magnetic field cre- 
ates an electric dipole moment, and several applications 
have been proposed for such responses. Such responses 
are observed in a variety of materials and from a variety 
of mechanisms. 1,2 From a theoretical point of view, the 
most intriguing part of the polarizability is that due to 
the orbital motion of electrons, because the orbital mo- 
tion couples to the vector potential rather than the more 
tangible magnetic field. 

The orbital magnetoelectric polarizability has also 
been studied recently in non-magnetic materials known 
as "topological insulators." These insulators have Bloch 
wavefunctions with unusual topological properties, that 
lead to a magnetoelectric response described by an E • 
B term in their effective electromagnetic Lagrangians, 3 
with a quantized coefficient. Qi, Hughes, and Zhang 3 
(QHZ) gave a formula for the coefficient of this term. 
For the specific case of topological band insulators, their 
result reproduces earlier formulas for the relevant topo- 
logical invariant, 15-17 but it is more generally valid: it 
describes a contribution to the magnetoelectric polariz- 
ability non just in topological insulators but in any band 



insulator. Their formula has a periodicity or ambigu- 
ity by e 2 /h that is related to the possibility of surface 
quantum Hall layers on a three-dimensional sample and 
generalizes the ambiguity of ordinary polarization. 

The same E • B coupling, known as "axion electro- 
dynamics" and originally studied in the 1980s, 4 was 
obtained in a previous paper by three of the present 
authors 5 using a semiclassical approach 6 to compute 
dP/dB, the polarization response to an applied magnetic 
field. However, in a general material, that semiclassical 
approach leads to an explicit formula for only part of the 
orbital magnetoelectric polarizability, the part found by 
QHZ. 3 The remainder, which is generically nonvanish- 
ing in materials that break inversion and time-reversal 
symmetries, is expressed only implicitly in terms of the 
modification of the Bloch wavefunctions by the magnetic 
field. 

In this paper, we develop a more microscopic approach 
that enables us to compute all terms in the orbital re- 
sponse explicitly in terms of the unperturbed wavefunc- 
tions, thereby opening the door to realistic calculations 
using modern band-structure methods (e.g., in the con- 
text of density- functional theory). Moreover, beyond its 
importance for computation, this expression clarifies the 
physical origins of the orbital magnetoelectric polariz- 
ability and resolves some issues that arose in previous 
efforts to describe the "toroidal moment" in periodic sys- 
tems. 

In the remainder of this introduction, we review some 
macroscopic features of the magnetoelectric response, 
while subsequent sections will be devoted mainly to a 
detailed treatment of microscopic features. The magne- 
toelectric tensor can be decomposed into trace and trace- 
less parts as 
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where a is traceless and 



is the trace part expressed in terms of the dimensionless 
parameter 8; a has the physical dimension of conduc- 
tance. The trace is the most difficult term to determine, 
as its physical effects are elusive. It should be noted that 
equality between dP l jdB^ and dMj/dEi only holds in 
the absence of dissipation and dispersion, which describes 
the low frequency, low temperature responses of an insu- 
lator. 7 ' 8 The placement of the indices in Eq. (1) is not 
essential for the arguments and calculations in this pa- 
per, and the reader can choose to treat a as a Cartesian 
tensor if desired. 9 As a Cartesian tensor, the traceless 
part decomposes further into symmetric and antisymmet- 
ric parts 

& ij = \ (<% + ) Otfj = \ ~ = -CijkTk, 
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where Tj = —eijkdijk/2 is the toroidal response. (Unless 
otherwise stated, in our work repeated indices are im- 
plicitly summed.) The terminology reflects that this part 
of the orbital magnetoelectric response is related to the 
"toroidal moment" , which is an order parameter that has 
recently been studied intensively; in a Landau effective 
free energy, the toroidal moment and the toroidal part of 
the magnetoelectric response are directly related. 10 ' 11 

The primary goal of this paper is to compute the con- 
tribution to a that arises solely from the motion of elec- 
trons due to their couplings to the electromagnetic po- 
tentials p4> and — j • A. We call this contribution the 
orbital magnetoelectric polarizability, or OMP for short. 
Other effects, such as those mediated by lattice distor- 
tions or the Zeeman coupling to the electron's spin, are 
calculable with known methods. 12 We shall only treat the 
polarization response to an applied magnetic field here; 
concurrent work by Malashevich, Souza, Coh, and one 
of us obtains an equivalent formula by developing meth- 
ods to compute the orbital magnetization induced by an 
electrical field. 13 

The magnetoelectric tensor's physical consequences 
arise through the bound current and charge, 3,8 given by 
Pb = — divP and Jb = d t P + curlM. Besides having 
a ground state value, each moment responds (instanta- 
neously and locally, as appropriate for the low-frequency 
response of an insulator) to applied electric and mag- 
netic fields, e.g., P % — Pq + Xe^3 + ol % aB^\ we will con- 
centrate on the magnetoelectric response. (Unless other- 
wise stated, in this article repeated indices are implicitly 
summed.) It is useful to allow a* , a material property, to 
vary in space and time by allowing the electronic Hamil- 
tonian to vary; this leads to a formula that covers the 
effects of boundaries and time-dependent shearing of the 
crystal, for example. Then the relevant terms are 

Jl = (aje <w - ay H )d k E x + (d t a))B> + e ijk {djo^Ei 
Pb = -afaBi - (d ia ))Bi, (4) 



We have used two of Maxwell's equations to simplify the 
first term in each line. The most important point to 
notice here is that ag does not appear except in deriva- 
tives, so that any uniform and static contribution to 6 
has no effect on electrodynamics. Hence in a uniform, 
static crystal, the components of a can be computed or 
measured from the current or charge response to spatially 
varying fields, given by the first term in each line. On 
the other hand, if we wish similarly to obtain ag from 
charge or current responses to applied fields, we need to 
consider a crystal that varies either spatially or tempo- 
rally, so that E or B will couple to diCtg or dtag, as in the 
second terms of Eqs. (4). These considerations, which we 
will elaborate later, motivate our theoretical approach to 
the OMP in this paper. 

We will proceed as follows. In Section II, we present 
the results of our calculation of the OMP in the 
independent-electron approximation. This section in- 
cludes a review of known results, followed by a discussion 
of the new contributions we compute and when those con- 
tributions can be expected to vanish (so that the OMP 
reduces to the form found in the literature previously). 
We follow these discussions with a detailed presentation 
of the calculation in Section III. This calculation involves 
a novel method for dealing with a uniform magnetic field 
in a crystal. An alternative derivation is presented in the 
Appendix. 



II. GENERAL FEATURES OF ORBITAL 
MAGNETOELECTRIC RESPONSE 

In this section we discuss properties of the OMP and 
its explicit expression in the independent electron ap- 
proximation. There is a natural decomposition into two 
parts, which is, however, not equivalent to the standard 
symmetry decomposition given in Eq. (I) of the Intro- 
duction. 

The first part is the scalar "Chern-Simons" term acs 
obtained by QHZ 3 that contributes only to the trace part 
ag. It is formally similar to the Berry-phase expression 
for polarization 14 in that it depends only on the wave- 
functions, not their energies, which explains the terminol- 
ogy "magneto-electric polarization" introduced by QHZ 
for acs- 3 The second part of the response is not simply 
scalar. It has a different mathematical form that is not 
built from the Berry connection, looking like a more typ- 
ical response function in that it involves cross-gap con- 
tributions and is not a "moment" determined by the un- 
perturbed wavefunctions. We label this term ac because 
of its connection with cross-gap contributions. This term 
does not seem to have been obtained previously although 
its physical origin is not complicated. 
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A. The OMP expression and the origin of the 
cross-gap term qg 

We first give the microscopic expression of the new 
term in the OMP and discuss its interpretation. The 
later parts of this section explain why the new term van- 
ishes in most of the models that have been introduced 
in the literature to study the OMP, and discuss to what 
extent the two terms in the OMP expression are physi- 
cally separate. The OMP expression that we discuss here 
will be derived later in Section III as follows: we com- 
pute the bulk current in the presence of a small, uniform 
magnetic field as the crystal Hamiltonian is varied adia- 
batically. The result is a total time derivative which can 
be integrated to obtain the magnetically induced bulk 



polarization. 

The most obvious property of the new term a G in the 
response is that, unlike the Chern- Simons piece, it has 
off-diagonal components; for instance, 8P x /dB y ^ 0. 
To motivate the expression for a G intuitively, we note 
that it is very similar to what one would expect based 
on simple response theory: An electric dipole moment, 
er, is induced when a magnetic field is applied. This field 
couples linearly to the magnetic dipole moment (e/4)(r x 
v — v x r) (this form takes care of the operator ordering 
when we go to operators on Bloch states). The expression 
we actually get for the OMP is expressed in terms of 
the periodic part of the Bloch wave functions u;k and 
the energies -E/k describing the electronic structure of a 
crystal: 
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Here the Berry connection A„ n ,(k) — i{u n k\dk a u n 'k) 
is a matrix on the space of occupied wave functions 
u nk , and the derivative with an upper index d a = dk a 
is a fc-derivative, as opposed to the spatial derivative 
di in Eq. (4). The velocity operator is related to the 
Bloch Hamiltonian, fiu*(k) = <9 4 _ff k , while the opera- 
tor / k is defined as the derivative d % V\s. of the projec- 
tion V onto the occupied bands at k. This operator is 
closely related to the position operator; its "cross-gap" 
matrix elements between occupied and unoccupied bands 
are (u mk |/ k Kik) = {u n k\ fl\u m u)* = -i(u m k\r z \u nk ), 
while its "interior" matrix elements between two occu- 
pied bands or two unoccupied bands vanish. Finally, the 
operator H' is introduced for generality, as discussed in 
Section III A; it vanishes for the continuum Schrodinger 
Hamiltonian and for tight-binding Hamiltonians whose 
hoppings are all rectilinear, and so will be ignored for 
most of the analysis that follows. Neglecting this sub- 
tlety, the form of a G is nearly what would be expected 
for the response in electric dipole moment to a field cou- 
pling linearly to the magnetic dipole moment. In the 
derivation presented in Section III, the term a G appears 
in abbreviated form at Eq. (37), and a G s follows imme- 
diately from Eq. (42). 

The main difference between the explicit form of a G 
and the naive expectation from the dipole moment ar- 
gument above is that a G excludes terms of the form 
(n|r|m)(m|v|n')(n'|r|n), for example, that include inte- 
rior matrix elements of r. In some sense, this omission 



is compensated for by the extra factor of 2 relative to 
the naive expectation and by a remainder term, namely, 
aesi the Chern-Simons part. The Chern-Simons term 
a G s alone has appeared previously. 3 ' 5 The next subsec- 
tion reviews the properties of a G s and gives a geometrical 
picture for its discrete ambiguity, which is not present in 
the a G term. We then explain how the existence of the 
previously unreported a G can be reconciled with previ- 
ous studies on model Hamiltonians that found only a G s, 
and then show that the two terms are more intimately 
related than they first appear. 



B. The Chern-Simons form, axion 
electrodynamics, and topological insulators 

The Chern-Simons response a G s has been discussed 
at some length in the literature. 3 ' 5 It does not emerge as 
clearly as a G from the intuitive argument above about 
dipole moment in a field; rather, in Ref. 5, it was derived 
by treating the vector potential as a background inhomo- 
geneity and utilizing a general formalism for computing 
the polarization in such a background. 6 

The most important feature of the microscopic expres- 
sion for the isotropic OMP is that it suffers from a dis- 
crete ambiguity. The dimensionless parameter 9 quanti- 
fying the isotropic susceptibility contains the term 
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which is only defined up to integer multiples of 2ir. This is 
tied to a "gauge" invariance: ground state properties of a 
band insulator should only be determined by the ground 
state density matrix p^, which is invariant under unitary 
transformations U nn >(\s.) that mix the occupied bands. 
Now, the Berry connection A is not invariant under such 
a transformation, but there is no inconsistency because, 
in the expression for 9c s, all the terms produced by the 
gauge transformation cancel except for a multiple of 2n. 
An analogous phenomenon, slightly easier to understand, 
is found in the case of electric polarization 14 
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which has invariance only up to a discrete "quantum," 
or ambiguity, which counts the number of times £/(k) 
winds around the Brillouin zone (e.g., if Un — e lk * a and 
Uu = 1, i ^ 1, then P x changes by one quantum). The 
Chern-Simons response acs behaves similarly, although 
the "winding" that leads to the ambiguity is more com- 
plicated (in particular, it is non- Abelian) . 

These ambiguities can be understood from general ar- 
guments, without relying on the explicit formulae. In the 
case of the polarization, the quantum of uncertainty of 
P x , e/S x , depends on the lattice structure, with S x the 
area of a surface unit cell normal to x. The ambiguity 
results because the bulk polarization does not completely 
determine the surface charge: isolated surface bands can 
be filled or emptied, changing the number of surface elec- 
trons per cell by an integer. For the magnetoelectric re- 
sponse, the quantum of magnetoelectric polarizability is 
connected with the fact that 9 gives a surface Hall con- 
ductance, as can be seen from the term J 5 = ( Vag) x E in 
Eq. (4). Therefore, the ambiguity in ag is just e 2 /h, the 
"quantum of Hall conductance," because it is possible to 
add a quantum Hall layer to the surface. (This remains 
a theoretical possibility even if no intrinsic quantum Hall 
materials have yet been found.) 

Now let us show that this ambiguity afflicts only the 
trace of the susceptibility. This can be seen directly by 
measuring the bound charge and currents. For example, 
all the components of a can be deduced from a measure- 
ment of pb in the presence of a nonuniform magnetic field 
[see Eq. (4)], but ag itself does not determine any bulk 
properties. 

More concretely, one can derive the ambiguities in the 
magnetoelectric response from the ambiguities in the sur- 
face polarization. In a periodic system, which for simplic- 
ity we take to have a cubic unit cell, the smallest magnetic 
field that can be applied without destroying the period- 
icity of the Schrddinger equation corresponds to one flux 
quantum per unit cell, or B = h/(eS), where S is again 
a transverse cell area. The ambiguity in the polarization 
of the system in this magnetic field corresponds to an 
ambiguity in dP/dB of 
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FIG. 1: A supercell admits a small magnetic flux, and the 
quantum of polarization transverse to the long direction is 
correspondingly small, but the quantum for polarization along 
the long direction is much larger. 



Hence on purely geometrical grounds there is a natural 
quantum e 2 /h of the diagonal magnetoelectric polariz- 
ability. 5 

In order to see that this uncertainty remains the same 
when a small magnetic field is applied (after all, a is 
defined as a linear response), we will have to construct 
large supercclls in a direction perpendicular to the ap- 
plied B (Fig. 1). While a supercell of N fundamental cells 
has a less precisely defined polarization (the quantum 
decreases by a factor N, so the uncertainty increases), 
the minimum field that can be applied also decreases by 
this factor, so that the uncertainty in the polarizability 
dP t /dB z (no sum) remains constant. On the other hand, 
if we consider the off-diagonal response, we can consider 
a supercell with its long axis parallel to the applied B. 
In this case, the polarization quantum remains constant 
as the supercell grows large and the minimum applied 
flux becomes small; the quantum in dP % JdB^ (for i ^ j) 
then becomes large, which means that the uncertainty 
vanishes. For this geometry, a small B acts like a contin- 
uous parameter, and the change in polarization induced 
by B can be continuously tracked, even if the absolute 
polarization remains ambiguous. 

Thus, with or without interactions, there is a funda- 
mental difference between the isotropic response and the 
other components of the response. For the trace-free 
components, we indeed do not find a quantum of uncer- 
tainty in the polarizability formula. In particular, if the 
toroidal response is defined by T, = —£^^(5^/2, then we 
believe that a "quantum of toroidal moment" 11 can only 
exist when there is a spin direction with conserved "up" 
and "down" densities. (This toroidal moment is typically 
defined as t = (1/2) Jrx fidr, with fi the magnetization 
density, 10 or more generally in terms of a tensor T lJ such 
that diT l i = — 2/i- 7 . 11 ) It then reduces to the polarization 
difference between up and down electrons. 

A particular class of materials for which the ambigu- 
ity in ag is extremely important is the strong topological 



insulators, 15-17 in which 9 = tt (Rcf. 3). These are time- 
reversal (T) symmetric band insulators. At first blush, T 
invariance should rule out magnetoelectric phenomena at 
linear order, since M and B are T-odd. However, the am- 
biguity by 2tt in 9 provides a loophole, since — n is equiv- 
alent to 7r. Here we regard the ambiguity/periodicity of 
9 as a consequence of its microscopic origin (alternately, 
its coupling to electrons); because 9 can be modified by 
2im by addition of surface integer quantum Hall layers, 
only 8 modulo 2tt is a meaningful bulk quantity for sys- 
tems with charge-e excitations. This is consistent with 
the gauge-dependence of the integral for acs- An alter- 
nate approach is to derive an ambiguity in 9 by assuming 
that the U(l) fields are derived from a non-Abelian gauge 
field. 4 The view here that periodicity of 9 results from the 
microscopic coupling to electrons is similar to the conven- 
tional understanding of the polarization quantum. 
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FIG. 2: Schematic band structure that leads to vanishing aa- 
The bands below the chemical potential are degenerate with 
energy E^_, while the bands above the chemical potential have 
energy = const. — 



C. Conditions causing aa to vanish 

It is worthwhile to understand in more detail the con- 
ditions under which the response ag is allowed. It is for- 
bidden in systems with inversion (P) or time-reversal (T) 
symmetry, which can be seen explicitly from the presence 
of three fc-derivatives acting on gauge-invariant matrices 
in the formula written in terms of Vk and i?k- 18 However, 
this alone is not sufficient to explain why ac did not ap- 
pear in the T-breaking models previously introduced to 
study the OMP. 3 - 5 - 19 This is explained by the fact that 
the intcrband contribution ac [Eq. (5b)] will also van- 
ish if dispersions satisfy the following "degeneracy" and 
"reflection" conditions: 

• At a given k, all the occupied valence bands have 
the same energy E£. 

• Similarly, all the unoccupied conduction bands 
have the same energy E^. 

• E£ + is independent of k (and can be taken to 
be zero). 

This can be seen immediately in an expanded form of the 
integrand of ac, [see Eqs. (Allc) and (Alld)] 
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where |n) = |w n k) and E n = E nk , etc. Such a structure is 
automatic when only two orbitals (with both spin states) 
are taken into account and the system has particle-hole 



and PT symmetries. PT symmetry guarantees that the 
bands remain spin-degenerate even if spin is not a good 
quantum number. To see this, recall that T acts on wave 
functions as io~ v K and maps k — > — k. Here, K is com- 
plex conjugation and a y takes the form of the usual Pauli 
matrix in the z basis of spin. Then P maps k — > k 
again, so that PT effectively acts as "T at each k." 20 
Then particle-hole symmetry implies that the dispersion 
is reflection-symmetric, E£ = — ££. 

Most model Hamiltonians discussed in the literature 
that access the topological insulator phase, 3 ' 5,15 ' 19,21 as 
well as the Dirac Hamiltonian (in the context of which 
the axion electrodynamics was first discussed 4 ), can be 
defined in terms of a Clifford algebra, 22 and this ensures 
that the dispersions arc degenerate and reflection sym- 
metric. The only exception of which we are aware is the 
model of Guo and Franz on the pyrochlore lattice, which 
has four orbitals per unit cell. 23 The topological insulator 
phase itself will not have a contribution from aa, since 
it is T- invariant, and so the Guo and Franz model will 
not show such a response; however, the addition of any 
T-breaking perturbation to their model should produce 
off-diagonal magnetoelectric responses. 

Finally, there is a simple mathematical condition that 
will cause aa to vanish. Namely, ag decreases as the gap 
becomes large without changing the wave functions, and 
in the limit of infinite bulk gap the only magnetoelectric 
response comes from the Chern-Simons part, which is not 
sensitive to the energies and depends only on the electron 
wave functions. 



D. Is the Chern-Simons contribution physically 
distinct? 



Apart from the ambiguity in acs that is not present in 
aa, there seems to be no real physical distinction between 
the two terms of the linear magnetoelectric response. We 
discuss two aspects that relate to this observation below. 

Localized vs. itinerant contributions 
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FIG. 3: A tetrahedral tight-binding molecule for spinless elec- 
trons, with one orbital per site and complex hoppings. The 
hopping integrals are all equal, except that those around one 
face have a phase of i relative to the other three. There are 
then two pairs of degenerate levels. 



The ambiguity in Oqs can be interpreted as a man- 
ifestation of the fact that bulk quantities cannot deter- 
mine the surface quantum Hall conductance, since a two- 
dimensional quantum Hall layer could appear on a sur- 
face independent of bulk properties. This suggests, per- 
haps, that the Chern-Simons term appears only in bulk 
systems with extended wave functions, and is a conse- 
quence of the itinerant electrons, while ac is a localized 
molecule-like contribution. However, this turns out not 
to be the case. 

Consider a periodic array of isolated molecules, which 
is an extreme limit of the class of crystalline insulators. 
Such a system has flat bands, with energies equal to the 
energies of the molecular states, since the electrons can- 
not propagate. It is certainly possible to construct a 
molecular system where all the unoccupied states have 
one energy and all the occupied states have another, by 
tuning the potentials. In this case ac will vanish. How- 
ever, such a molecule can still display a magnetoelectric 
response; it will therefore have to be given by acs (and so 
restricted to diagonal responses). For example, consider 
the "molecule" of Fig. (3) with the shape of a regular 
tetrahedron. If the two low-energy levels are occupied, 
the magnetoelectric response is 



±51 



V6 fi ' 



where P l here is the electric dipole moment divided by 
the volume of the tetrahedron; the sign of the polariz- 
ability reverses when the complex phases are reversed. 
This shows that the Chern-Simons term does not require 
delocalized orbitals. 
Additivity 

Another argument against distinguishing between the 
Chern-Simons part and the rest of the susceptibility is 
based on band additivity. When interactions are not 



taken into account, each occupied band can be regarded 
as an independent physical system (at least if there are no 
band crossings). Applying a magnetic field causes each 
band n to become polarized by a certain amount P n , and 
so the net polarization should be P = J2 n The Pauli 
exclusion principle does not lead to any "interactions" be- 
tween pairs of bands, because the polarization (like any 
single-body operator) can be written as the sum of the 
mean polarization in each of the orthonormal occupied 
states. 

Now the Chern-Simons form does not look particu- 
larly additive in this sense, and is not by itself. Be- 
cause it is the trace of a matrix product in the occu- 
pied subspace, it necessarily involves matrix elements be- 
tween different occupied states, while an additive formula 
would not. Nevertheless, ac and acs are together ad- 
ditive, as can be seen most simply in Eq. (A9), where 
the two terms combine into a single sum over occu- 
pied bands. In terms of qg and acs separately, one 
finds that when the values of ac, assuming just band 
1 or 2 is occupied, are added together, some terms oc- 
cur that are not present in the expression for ac{\ + 2) 
(where both bands are occupied), and vice- versa. Us- 
ing Eqs. (Allc) and (Alld) we see, in fact, that ac is 
a sum of contributions which depend on three bands, as 

«G =E„,m,m' C ( n ; TO > TO ')+E„,„',m Z) ( n >™'; m )- TermS 

such as C(l; 2,m') are not present in the expression for 
«g(1 + 2). (Likewise D(l, 2;m) appears in ac(l + 2) 
but not in ac(l) and «g(2).) Adding up the discrepan- 
cies, one finds that the energy-denominators all cancel, 
and the non-diagonal terms from the Chern-Simons form 
appear! 

Seemingly paradoxical is the fact that for band struc- 
tures satisfying the degeneracy and reflection conditions 
of the last subsection, the magnetoelectric susceptibility 
is given by the Chern-Simons term alone, which does not 
seem to be additive. However, the additivity property ap- 
plies only to bands that do not cross. It does not make 
any sense to ask whether the susceptibility is the sum 
over the susceptibilities for the systems in which just one 
of the degenerate bands is occupied, since those systems 
are not gapped. 



III. THE OMP AS CURRENTS IN RESPONSE 
TO CHEMICAL CHANGES 

Now we will tackle the problem of deriving the formula 
for the OMP a discussed in the last section. There are 
two impediments we need to overcome, a physical one, 
and a more technical one (which we will overcome start- 
ing from an insight of Levinson). 24 

In order to determine a, we would like to carry out a 
thought experiment in which a crystal is exposed to ap- 
propriate electromagnetic fields. For specificity, we will 
apply a uniform magnetic field. To make the calculation 
of the response clean, we wish to deal with an infinite 
crystal. Then the polarization does not simply reduce 
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FIG. 4: As outlined in the text, (a) turning on a magnetic 
field produces a macroscopic polarization through the flow of 
surface currents, while (b) varying the crystal Hamiltonian in 
the presence of a fixed magnetic field produces a polarization 
through the flow of current through the bulk. 



to the first moment of the charge density, 14 so we will 
instead have to calculate the current or charge distribu- 
tion induced by the fields, and then use Eq. (4) to deduce 
a. If both the crystal and the electromagnetic fields are 
independent of space and time, there is no macroscopic 
charge or current density. We will assume spatial unifor- 
mity, so that there are two choices for how to proceed. 
Either the magnetic field can be varied in time or the 
crystal parameters, and thus a, can be varied. In either 
case, we measure the current that flows through the bulk 
and try to determine a. As ever, the diagonal response 
olq is the most difficult to capture: while either time- 
dependent experiment can be used to determine a, only 
the latter approach sheds light on the value of ag . 

To see why otg can be determined only in this way 
(given that we want to work with a spatially homoge- 
neous geometry) , let us discuss how currents flow through 
the crystal. The necessity of varying the crystal in time 
can be deduced from Maxwell's equations (see below) 
but we will give a more intuitive discussion here. Sup- 
pose that 5 = 0. Then in an applied magnetic field 
there is a polarization P = agB: thus the crystal gets 
charged at the surface. As the magnetic field is turned 
on, this surface charge has to build up (charge density 
n • P). This occurs entirely due to flows of charge along 
the surface. Suppose, for example, that the sample is a 
cylinder (radius R) with the magnetic field along its z 
axis, as illustrated in Fig. 4(a). Then an electric field 
Emd = —BRcp/2 is induced at the surface according 
to Faraday's law. Besides being the magnetoelectric re- 
sponse, 9 also represents the Hall coefficient for surface 
currents. Therefore, a current of J s = otgBR/2 flows to 
the top of the cylinder, adding up to a surface charge of 
2ttR J J s (t)dt — agBfirR 2 and producing the entire po- 



larization agBf. No current flows through the bulk! In 
fact, the Hall conductance on the circular face produces 
a radial current as well, so that the charge distributes 
over the surface rather than just accumulating in a ring. 
Note that the surface current grows with the radius of the 
cylinder. This sounds like a nonlocal response, but it can 
be understood as follows: the electric field is determined 
by the non-local Faraday law, but the crystal's response 
to the electric field (namely, the surface current) is local. 

The current distribution can be understood directly 
from Maxwell's equations: there are two contributions 
to the bulk current, dtP and curl M. The polariza- 
tion is ag~B while the magnetization is indirectly pro- 
duced by the induced electric field, agE^d- The two 
contributions thus cancel by Faraday's law in the bulk: 
jbuik _ Q, e (g t B_)_ cur l E) = 0. There is a surface current 
because ag is discontinuous there. 

On the other hand, if changes in time, while the mag- 
netic field is time- independent (as in Fig. 4b), the polar- 
ization at the ends of the cylinder builds up entirely by 
means of flows of charge through the bulk. Surface flows 
cannot be large enough to explain the net polarization 
in this situation. Since there is no induced electric field, 
the surface current is just proportional to the lateral sur- 
face area and is negligible compared to the bulk current. 
Therefore the bulk current is equal to (dtag)B and can 
be integrated to give agB. 

For the other component of the OMP, a, cither thought 
experiment can be used. The simplest approach, how- 
ever, is still the crystal-variation method, since the sur- 
face currents are negligible in that case, 25 and in any 
case this method allows us to find all the components of 
a simultaneously. 

Difficulties with the operator r and uniform magnetic 
fields 

There are two technical difficulties in the theory. First, 
the operator r has unbounded matrix elements and thus 
the matrix elements of the magnetic dipole moment 
(e/4)(v x r — r x v) are not well-defined. This rules out 
the straightforward use of perturbation theory to calcu- 
late the electric dipole moment of an infinite crystal in a 
uniform magnetic field. Second, if we consider a crystal 
in a uniform magnetic field, Bloch's theorem does not 
hold. Although the magnetic field is uniform, the vector 
potential that appears in the Hamiltonian depends on r. 

We avoid the problems of r as follows. The key idea 
is to work with the ground state density matrix, rather 
than wave functions. The individual eigenstates change 
drastically when a magnetic field, no matter how small, 
is applied (consider the difference between a plane wave 
and a localized Landau level). However, the density ma- 
trix of an insulator summed over the occupied bands only 
changes by a small amount when B is applied; over short 
distances the magnetic field cannot have a strong effect 
(even in the example of Landau levels) , and the density 
matrix has only short-range correlations because it de- 
scribes an insulating state. More technically, we show 
(subsection III A) that the broken translation invariance 
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of any single-body operator O (such as the density ma- 
trix) can be dealt with by factoring out an Aharonov- 
Bohm-like phase from its matrix O rir2 . This solves the 
problem of the nonuniform gauge field and leads to ex- 
pressions that depend only on differences between r's. In 
addition, since the exponentially decaying ground-state 
density matrix appears multiplying every expression, the 
factors of ri — r 2 arc suppressed. 

The calculation then proceeds as follows. First, using 
the symmetries of the electron Hamiltonian in a uniform 
magnetic field, we find how the density matrix changes 
in a weak magnetic field. Next we compute the current 
response to an adiabatic variation of the crystal Hamilto- 
nian. Finally, we show that this current can be expressed 
as a total time derivative, and therefore can be integrated 
to give the polarization; at linear order in B we can read 
off the coefficients, the magnetoelectric tensor a. 



This argument shows how to couple general Hamiltoni- 
ans to uniform fields: H = exp[(«e//i) d£- A][H 0ri r 2 + 
H' rir2 (B)]. The vector potential appears explicitly only 
in A, while H'{B) gives the rest of the dependence on 
the magnetic field. The Schrodinger Hamiltonian (10) is 
obtained if we take 



H, 



0,rir 2 



5<*Hv 2 -r,). (15) 



and set H 1 = 0. Our results also apply to tight-binding 
models. We introduce H' to capture the possibility that 
in a tight-binding model the hoppings will not be rec- 
tilinear, and hence that the phases in Eq. (14) do not 
capture the full field dependence of the Hamiltonian. 



B. The ground state density operator 



A. Single-body operators for a uniform magnetic 
field 



Recall the form of the Schrodinger Hamiltonian for a 
single electron in a crystal and under the influence of a 
magnetic field, 

i/s(p ' r) = i [p ~ eA(r)]2 + y(r) ' (10) 

where V(r + R) = V(r) for lattice vectors R. The neces- 
sity of using the vector potential A seems at first to spoil 
the lattice translation symmetry one would expect in a 
uniform magnetic field. However, as noted by Brown 26 
and Zak, 27 a more subtle form of translation symmetry 
remains. In particular, choosing the gauge 



A=-Bxr, 
2 



(11) 



the Hamiltonian has "magnetic translation symmetry" : 

H s (p, r + R) = e ieB < Rxr ^ 2H H s {p, r ) e -«B.(Rxr)/2K_ 

(12) 

This condition defines magnetic translation symmetry 
for general single-body operators. Any operator O pos- 
sessing this symmetry can be written in the position basis 
as 



/n _ /n -ieB-(riXr 2 )/2ft 
v - / rir 2 — v - / rir 2 t ; i 

where O has lattice translation invariance, 



O 



ri+R,r 2 +R — 



(13a) 



(13b) 



Note that the phase is just (ie/h) J d£-A calculated along 
the straight line from r 2 to ri, which agrees with the 
intuition that comes from writing the second-quantized 
form of the operator, 



0= I d 3 r 1 d 3 r 2 O rir2 cl i e 



t -ieB-(nxr 2 )/2fi 



(14) 



We find it convenient to work with the one-body den- 
sity matrix p 9 , or equivalently the projector onto the oc- 
cupied states, whenever possible, because it is a basis- 
independent object. Also, in an insulator, p 9 ir2 is expo- 
nentially suppressed in the distance |r 2 — ri|, which tem- 
pers the divergences that arise from the unboundedness 
of r. 28 In any case, if the ground state is translationally 
symmetric, the structure described above will apply to 
p g and we can be sure that the density matrix has trans- 
lational symmetry apart from a phase: 



P 9 



-»eB-(ri xr 2 )/2ft 



(16) 



where p 9 possesses the translation symmetry of the crys- 
tal lattice and hence should connect smoothly to the field- 
free density matrix. Hence we will write 



p 9 =p + p', 



(17) 



where po is the density operator of the crystal in the 
absence of the magnetic field. 

Density matrix perturbation theory: Now we have to 
calculate p', using a kind of perturbation theory that 
focuses on density matrices rather than wave functions, 
since the wave-functions suffer from the problems dis- 
cussed above. This perturbation theory starts from two 
characteristic properties of the density matrix: it com- 
mutes with H, and for fermions at zero temperature it is 
a projection operator. The latter means that all states 
are either occupied or unoccupied, so the eigenvalues of 
the density operator are and 1, which is formalized as 



p 9 p 9 = p 9 



(18) 



(idempotency). 9 Expressed in the position basis, 



9 = [ drop 9 p 9 

9 - rfr^n 9 7,3 p -(*e/2fi)B-(n Xr 2 +r 2 xr 3 +r 3 xn) 

rir 3 / ul 2Frir 2 Fr 2 r 3 c 



(19) 
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The exponent is just —i(t>i23/4>o, proportional to the mag- 
netic flux through triangle 123, and the exponential can 
be expanded for small B. At first order this gives 

Pr ir3 = J dr2[/>r ir2 P0r 2 r 3 + P0r ir2 Pr 2 r 3 

P0rir 2 P0r 2 r 3 (#123/^0)]. (20) 

The problem of the unbounded r's is resolved in this equa- 
tion because the area A123 of the triangle is finite and 
independent of the origin, and also suppressed by the 
factor of p. 

Calculation of p': In the last term of Eq. (20), we can 
rewrite 2^123 = ri x r 2 + r 2 x r 3 + r 3 x ri = (r 2 - ri) x 
(r 3 - r 2 ) and then use (r 2 - ri)ponr 2 = [po,r]nr 2 , etc., 
to obtain 

(1 -Po)p'(l -Po) - Pap'po = - { ^ B - dPo,r] x [Po,r]). 

(21) 



If we define 



H = H Q + H', 



(22) 



then analogous manipulations (including (ri— r2)H rir2 = 
ifiv rit! ) on the equation [H, p 9 ] = give 

[p',H ] = • ([p ,r] x v - v x [p ,r]) - [p 0} H'}. (23) 

Eqs. (21) and (23) have an intuitive meaning. The for- 
mer equation determines the "interior" matrix elements 
of p' , those between two occupied or two unoccupied 
states of the zero-field Hamiltonian. An perturbation 
with the full crystal symmetry does not change the in- 
terior matrix elements of the density matrix because of 
the exclusion principle. 30 In our case, however, multiply- 
ing po by the phase e ( ie / 2fi ) B ' r i xr 2 gives a density matrix 
with the correct magnetic translation symmetry, but also 
changes the momentum of the states and so results in a 
small probability for states to be douby occupied. There- 
fore p' must correct for this "violation of the exclusion 
principle." On the other hand Eq. (23) determines the 
"cross-gap" matrix elements of p' (those between unoc- 
cupied and occupied states). These matrix elements cap- 
ture the expected "transitions across the gap" induced 
by the field. The rest of this section is devoted to calcu- 
lating all these matrix elements. The results are given in 
Eqs. (24) and (28); the derivations could be skipped on 
a first reading. 

Calculation of p' . Precisely speaking, Eq. (21) gives 
the matrix elements of p 9 between pairs of occupied (n 
and n') or unoccupied (m and m') states: 



<V>mk|(l - p 9 )\An^) = 6 mm , - —B>e jab F% m ,(k), 



Ah 



where T is the non-Abelian Berry curvature associated 
with the occupied bands, 

Tt n , = i(u nk \d a V u d b Vu - d b V k d a Vu\u n ,u) 

= d a A b nn ,-d b A a nn ,-i[A a ,A b } nn ,, (25) 

and JF is the corresponding quantity for the unoccupied 
bands. To derive these relations, we use 



BZ 



d 3 k 

(2tt) 3 



"Pke 



(26) 



where V — 5^ nocc \ u nk)( u nk\ is the projector onto filled 
bands at k. This gives 



-/ 

Jbz 



d 3 k 

BZ (27T) 3 

d 3 k 



*[pV] = / 7^T3 elk ""( V k^ k )e- lkr 



(27) 



after discarding a total derivative. The notation / = 
VkP was introduced in Eq. (5). 

By contrast, Eq. (23) describes to what extent p 9 fails 
to commute with H n , the crystal Hamiltonian, and gives 
the matrix elements of p' between occupied and unoc- 
cupied states. In this sense it is analogous to the more 
usual results for density-matrix perturbation theory. 30 In 
the basis of unperturbed energy eigenstates, 



, , , ,, , v . e • Kk|{d a Pk,d fc g k }K k ) 

{lpnk\P IV'mk) = l^B 3 e jab p p 



(28) 



(24) 



Recall that hv = d h H-^ and that H' is introduced only to 
capture unusual situations such as tight-binding models 
with non-straight hoppings, and vanishes for the contin- 
uum problem. Eqs. (24) and (28) are the key technical 
results of this formalism, good to linear order in the mag- 
netic field. 



C. Adiabatic current 

Now we need to calculate the current as the Hamilto- 
nian is changed slowly as a function of time, as in the 
ordinary theory of polarization. 14 ' 31 We have to be care- 
ful, however, since the current vanishes in the zero-order 
adiabatic ground state described by density matrix p 9 (t). 
It is necessary to go to first order in adiabatic perturba- 
tion theory, which takes account of the fact that the true 
dynamical density matrix p(t) has an extra contribution 
proportional to dH/dt = H. However, once the current 
has been expressed in terms of p, which is proportional 
to H, the distinction is no longer important and the adi- 
abatic approximation can be made. 
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Preparing for the adiabatic approximation: We can 
write the current as 

3(t) = ^Trp(t)v=^ l -Trp[H,r] (29) 

where £1 is the crystal volume. Here H is the full 
Hamiltonian including the magnetic field. By unitar- 
ity of time evolution, it remains a projector if the ini- 
tial state describes filled bands only. The operator r ap- 
pears here, but in a commutator. Since (ri | [O, r] |r 2 ) = 
(r 2 — ri)O rir2 , such expressions do not suffer from the 
difficulties of an "unprotected" r, namely its unboundcd- 
ness. We can only use cyclicity of the trace to the extent 
that this property can be preserved. In particular, the 
expression Trr[p, H], which seems formally equivalent to 
Eq. (29), poses problems, but 



Integrating the results: The four terms in the current 
can be collected and rearranged into the form 



J = Jg + Jcsi + JCS2 



(35a) 



and integrated with respect to time as follows. The first 
term in Eq. (34) can be rewritten with por 1 r 2 Por 2 r 3 p' r3ri = 

r 2 P0r 2 r 3 Pr 3ri ) A)rir 2 P0r 2 r 3 Pr 3 r ± P0r 1 r 2 p0r 2 r 3 Pr 3 r 1 

and combined with the next two terms to give 

J G = ^Tr[p ,r][p>o] (35b) 
Jcsi = -3-Trp'[p , [po,r]], (35c) 

while the final term in Eq. (34) (i.e., the term involving 
^123) becomes 



e i , 



JM = ^Tr[p, [p,r}}[p,H] 



(30) 



does not. This expression can be derived from Eq. (29) 
using again the idempotency of p (pp = p). Using the 
equation of motion for the density matrix, 



ihp{t) = [H(t), P (t)} 



(31) 



and making the approximation p m p 9 on the right-hand 
side at this stage, the current becomes 

J = h j dT ^ dr ^ - 2l "2 + r 3K ir2 /^ r3 P? 3ri . (32) 

Magnetic field dependence of the current: The consid- 
erations given in the last subsection make the integrand 



Pv 1 r 2 Pr 2 r 3 Pr 3 r 1 P 



- n 9 



rir 2 rr 2 r 3 rr 3 i"i v 



(33) 



where, again, </>i 23 = B- (ri x r 2 +r 2 x r 3 +r 3 x ri)/2 is the 
magnetic flux through the triangle with vertices rir 2 r 3 
and does not suffer from the pathologies of r itself, which 
allows us to expand e _4 ^iW0o i _ i0 123 /0 o to lowest 
order in B (recall again that the matrix elements of p arc 
exponentially suppressed with distances). 

Recalling our division p 9 = p + p 1 where p' is of first 
order in the magnetic field, Eq. (32) becomes 



A3rir 2 P0r 2 r 3 P r3ri 



Pr 1 r 2 P0r 2 T 3 P0r 3 T 1 + P0rir 2 Pr 2 r 3 /'0r3ri 



^123 



"P0rir 2 P0r2r 3 P0r3ri 



(34) 



at first order. The rest of the calculation involves substi- 
tuting the expressions for the magnetic-field dependence 
of p 9 obtained earlier, and integrating the result to ob- 
tain a. The energy-dependent part of a, namely oiq, will 
come from the mixing of the occupied and unoccupied 
bands, Eq. (28). The Chern-Simons term will come from 
the "exclusion-principle-correcting" terms, Eq. (24), as 
well as the </>i 23 term in the previous equation. 



3cS2 =- i xM Bie j"bTr[p ,r}[p ,r a }[r b ,po}+c.c. (35d) 

upon rewriting 

(ri - 2r 2 + r 3 )(ri x r 2 + r 2 x r 3 + r 3 x ri) 
= (ri - r 2 )[(ri - r 3 ) x (r 2 - r 3 )] 

+ (i"3 - i"2)[(ri - r 2 ) x (n - r 3 )]. 

The total derivative term Jg [Eq. (35b)] can be written 

J G = dt(a G ))B* (36) 

with ao as given in Eq. (5b), 

e 2 „ ^ (n\d i V\m)(m\e jab {d a H,d b V}\n) 



n occ 
m unocc 



E n — E r , 



+ 2eIm (Md^HimldH'/dB^n) ^ 



n occ 
m unocc 



E n — E n 



where the BZ integral and the dependence on k have 
been suppressed, and \n) — |u„), etc. This result fol- 
lows immediately upon taking the trace in the basis of 
energy eigenstates. Matrix elements of [po> rJ ] appear as 
= d l Vk, from Eq. (27), and the cross-gap matrix ele- 
ments of p' in are given in Eq. (28). Note that since Jg 
is a total time derivative, aa is uniquely defined for a 
given Hamiltonian (this assumes the existence of a refer- 
ence Hamiltonian with ciq = 0, that is, the existence of 
a topologically trivial, time-reversal-invariant band insu- 
lator). 

In J G S2 [Eq. (35d)], we can replace r ->• [[r, po], po] 
in the third commutator. This has the same cross-gap 
matrix elements as r; the interior matrix elements do not 
contribute to the trace because the other three factors, 
po and two components of [po,^], have only cross-gap 
matrix elements. Then 



Jc52 = i-^B 3 e jab Tr[po,r][p Q ,r a ][[[p ,r ], p ], po] + c.c. 
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or 

J C S2 = ~B^e jab TrV k r k d a r k [[d b r k ,P k ],r k } +c.c, 

where an integral over k is suppressed for brevity and 
the trace is taken in the Hilbert space at k. Dropping 
the subscripts k everywhere, this can be expanded and 
rearranged to give 

Jcs2 - ^e jab TrV{[vv,v}d a vd b v 

+ 2[V,d b V][VV,d a V] 
+ 3(Vd h Vd a VVV + d b VVVVd a V)}. (35d') 

In manipulating these strings of projection operators and 
their derivatives, it is very useful to realize that deriva- 
tives of projectors only have cross-gap matrix elements: 
Vd a VV = Qd a VQ = 0, where Q = l-V is the projector 
onto unoccupied bands. This means, for example, that 
V(VV)(V) = V(VV)Q(V)V. 

To Jcs2 we must add Jcsi [Eq. (35c)], 

Jcsi = ^rB j e jab Tr(V - Q)[P, VV]d a Vd b V 

= : ^-B^ jab Tir{[v,vr}d a rd b r 

- (Vd b Vd a VVr + d b VPVVd a V)}, (35c') 

to get 

Jcs = Jcsi + Jcs2 

p 2 

= —B j e jab TvV{[V,VV]d a Vd b V 

+ [P,d b V][W,d a V}}. (38) 

By checking the different components explicitly one can 
see that this is 

e 2 

Jcs - H—TrV{[P, d x V\ [d v V, d z V\ 

+ [P, d v T] [d z V, d x T] + [V, d z T] [d x V, d v V]}, (39) 
so we get the "topological current" 

J cs = -B- / — -tr(T tx T vz + P V T ZX + T tz T xy ), 
n Jbz 

(40) 

where the lower-case trace (tr) is only over the occupied 
bands, and the Brillouin-zone integral has been restored. 

It remains only to show that Jcs is a total time deriva- 
tive that integrates to acsB. Allowing the indices to run 
over t, x, y, z, in that order (so that e tX yz = +1), 

Jcs = -B— f ,r. e abc dtrJ rab J rcd 
8H J BZ (2tt) 3 

- - B f^' L (0 a "' r ( A "' fA * - i A>A ' Al ) ■ 

(41) 



The derivatives with respect to k x ,k y , k z will vanish when 
integrated over the Brillouin zone assuming that A is 
defined smoothly and periodically over the zone, leaving 
just 

Jcs = -B|U - *|^«) , 

. (42) 

where the indices now only run over xyz, as originally. 
This obviously gives acs as in Eq. (5c), completing the 
proof. It must be reiterated that this integral is not al- 
ways entirely trivial. In particular, if the adiabatic evo- 
lution brings the crystal back to its initial Hamiltonian 
in a nontrivial way, the Brillouin zone integral need not 
return to its initial value because A is not uniquely de- 
fined. In other words, J dtJcs can be multivalued as 
a function of the Hamiltonian deformation parameters. 
However, the change can only be such that 8 changes by 
an integer multiple of 2n, as discussed in subsection II B. 



IV. SUMMARY 

The theoretical calculation of the magnetoelectric po- 
larizability in insulators presents a difficulty similar to 
that known well from the theory of polarization; both 
quantities suffer an inherent ambiguity in the bulk. The 
magnetoelectric polarizability adds another level of dif- 
ficulty because the vector potential is unbounded and 
breaks lattice translation symmetry. However, we have 
developed a formalism that allows us to deal directly with 
a uniform magnetic field. In the appendix, we further 
show that a long-wavelength regularization of the vec- 
tor potential together with a suitable generalization of 
the polarization (to deal with the broken crystal symme- 
try) provides a (relatively) simple, though less rigorous, 
way to compute the response function. The final expres- 
sion for the OMP rederives known results for particular 
model systems and topological insulators and completes 
the picture with additional terms that have a relatively 
straightforward and intuitive interpretation. We hope 
that these results and the method of their derivation will 
be valuable for future work on magnetoelectric effects 
and topological electronic phases. 

The authors gratefully acknowledge useful discussions 
with S. Coh, A. Malashevich and I. Souza. The work 
was supported by the Western Institute of Nanoelectron- 
ics (AME), DARPA OLE (AMT), NSF DMR-0804413 
(JEM), and NSF DMR-0549198 (DV). 

Appendix A: Calculating the OMP using Static 
Polarization 

As noted in the text, matrix elements of the opera- 
tor r are ill-behaved in a basis of extended, Bloch-likc 
states. That problem was solved by working with the 
density operator p, whose matrix elements are exponen- 
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tially suppressed with distance. Another approach is to 
use a Wannier-like basis of localized states. In this ap- 
pendix, we take this approach to present an alternative 
derivation of the OMP. 

The Bloch functions VVik(r) of the unperturbed crystal 
will evolve, under the application of a long-wavelength 
magnetic field A = A sin q • r, into the exact energy 
eigenfunctions ^„k(r). These no longer have a sharp 
crystal momentum k, but may be expanded in a per- 
turbation series in the unperturbed VVik( r )- Then the 
analogue to the standard Wannier function w nR (r) for 
lattice vector R will be 



W nR (r) = 4= 



d 3 k 



*nk(r)e 



-ik R 



= w nR (r) + 5w nR (r), 



(Al) 



where f2 is the volume of the crystal and N is the number 
of unit cells. The Wannier orbitals centered at R become 
polarized when the magnetic field is applied, and this 
distortion gives a polarization density of 



5P(R) = ( w nR\er\Sw nR ) + c.c. 



(A2) 



Although it is not obvious that the bulk polarization ap- 
pearing in Maxwell's equations is the same as the po- 
larization of a set of Wannier orbitals, this expression 
leads to Eq. (5). To ensure that the Wannier orbitals 
are localized, we will have to suppose that each band has 
a vanishing Chern number, 32 so that the phase of w„k 
can be chosen so that it is a periodic function of k. In 
this case the unperturbed Wannier functions are local- 
ized, and (though there are usually subtleties in defin- 
ing Wannier functions in a magnetic fields), 33 the reg- 
ularization used here leads to localized orbitals. Pre- 
sumably these arguments can be extended to the case 
where the total Chern number for all occupied bands 
C ij = E„occ Jbz d 3 k^ n (k)/2n vanishes. 

Here we want to take a relatively direct approach to 
perturbation theory in the field, and write 34 



*«k = "0nk 



_ eB v 

~ 2^4" 



V'Zk+q 



(A3) 
(«->-«) • 



For definitcness we take A = — (B/q)s'm(qy)x, and 
the velocity operator can be alternatively expressed as 
v x = d x H^/h, with the Bloch Hamiltonian of the 
unperturbed crystal. 

Then the first-order correction to the dipole moment 
of the generalized Wannier functions will be 

fF*(R) = e ]T / dr I (Vv k -( R -)< k (r)) 

x / 7^^nk'(r)e- 4k '- R + c.c. (A4) 



The position integral must be taken over the whole crys- 
tal at this point. In the integral over k, r l can be con- 
verted into a k derivative of the exponential, and then 
partial integration leaves a factor —id^u* (the boundary 
term vanishes because the Bloch function ip is strictly 
periodic in k). Then 



5P l = - 



e 2 B 
2q 



E 

n occ 
I 



(d l U„u | M ik ) (UIU I V X I U„k-q) 



E. 



rak- 



Eik + it 



k-q/ gl q.R 



-((/->• -9) 



+ c.c. (A5) 



(From now on, we will omit the integral over k and 
the associated factor of (27r) 3 .) Because of the varia- 
tion of the magnetic field the magnetoelectric polariza- 
tion 6P % (R) — aJS J '(R) should vary as cosq • R. The 
polarization seems to have both cosine and sine terms, 
but the coefficient of the latter is — B sin(qy)C lx / q, and 
the vanishing of C is a prerequisite for using Wannier 
functions. 

To lowest order in q and B, then, the magnetoelectric 
response is 

i _ _ e 2 Q (l9'n n k|MZk)(^kK|Unk-q) 

«i- 2 e ^2. E n ^-E lk + it +C ' C " 

n occ 

(A6) 

where we have symmetrized over Landau gauges to make 
the expression nicer. 

Switching to the shorthand \n) — |u„k), 



(d l n\l)(l\d a H\d h n) 
E n -Ei+ie 

(d l n\l)(l\d a H\n)d b E n 



(E n -Eh + itf 



(A7) 



Simplifying the second term of this expression makes use 
of the "Sternheimer equation" 

(8 a H)\n) = (E n - H)\8 a n) + (d a E n )\n) (A8) 

and the antisymmetry in the indices a and b to give 

1 {d i n\l){l\d a (H + E n )\ 8 b n) 



2h t3ab E 



n occ 
I 



E n -Ei+ ie 



+ c.c. (A9) 



Note the formal similarity to the expression for orbital 
magnetization, 

Mj = hmJ2e jab (d a n\(H + E n )\8 b n) 

n occ 

= -^E^^' 9 "^ + E n )\d b n), (A10) 



in particular the appearance of the combination H + 

-E„.. 34 ' 35 
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To bring our compact expression into the form given 
in terms of ag and acs m the main text, we need to 
break the sum over I into contributions from occupied 
and unoccupied states. Omitting the factor (e 2 /2H)ej a b 
for the moment, the sum over the occupied states takes 
the form 



(d*n\l){l\d a (H + E n )\d b n) 

2 ; " " ■ ■ +C.C. 



n,l 
occ 



E n -E[+ ie 



^ , (l\d a (H - E^n) + (d»l\d a (H - E n )\n) 
~^ dn W E n -E l+ ze 



n,l 
occ 



J2{d i n\n')(d b n'\d a n) 



(Alia) 



n.n 

occ 



using the antisymmetry in a and b and the Sternheimer 
equation again. Because the two sums are not symmetric 
when we take I in the unoccupied space, however, the 
terms do not cancel as nicely. Inserting a resolution of 
the identity, broken into two parts, gives: 

^ (d i n\m)(m\d a H\n')(n'\d b n) 

E J ri +c - c - 

{d t n\m){m\d a n l ){n'\d h n) + c.c. (Allb) 



n,n occ 
m unocc 



n,n occ 
m unocc 



^ (d*n\m)(m\d a n>)(n'\d b H\n) | ^ .^^ 



n,n occ 
m unocc 



E n - E n 



^ (d t n\m)(m\d a n)d b E n 
> - h C.C 



n occ 
m unocc 



E — E 



and 



E — ^ — 5 — +c - c - ( Alld ) 



n occ 
m,m' unocc 



F — F 



+ E 

n occ 
m unocc 



(d i n\m)(m\d b n)d a E n 
E — E 



+ C.C 



The unnumbered pieces of these equations cancel by an- 
tisymmetry in a and b. 

Defining V as the projector onto occupied bands as in 
the text, Eqs. (Allc) and (Alld) combine to give 



e 2 ^ (n\d*T\m)(m\{d a H,d b T}\n) 
[a G )j = -rz tjab 2^ S S + c - c -> 



E n - En 



m unocc 



(A12) 

which is equivalent to Eq. (5b) upon identifying v a with 
d a H and /P with d % V . This quantity has the crucial 
property that it is "gauge invariant," meaning that it 
can be written as a matrix trace, and hence does not 
change under a change of basis of the Hilbert space. Of 
course, this property is not evident here, where the for- 
mula makes explicit reference to energy eigenfunctions 
and their energies, but it follows from the expression in 
terms of a matrix given in Eq. (35b). The remainder, 
Eqs. (Alia) and (Allb), becomes 



(«cs)'j 



A a d b A c 



2i 



A a A b A c 



, (A13) 



which reproduces Eq. (5c). 
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